function h = plot_usef(grid,base,NU,xi,xie)
% h = plot_usef(grid,base,NU,xi,xie)
%
% Plot the u^* reconstruction, including the exponential factor.
% Assuming that NU contains the coefficients of nu with respect to
% base, u^* is given as
%  u = nu*exp(xi)
%
% xi: function pointer returning the value xi
% xie: vector with the correction values for xi on each element:
%       xi_eff|_K = xi - xie(K)
 
 meshg = load('plot_mesh.mat');

 % evaluate the basis on the meshg.pg nodes
 PHI = zeros(base.pk,size(meshg.pg,2));
 for i = 1:base.pk
   p_si = base.p_s{i};
   for j = 1:size(p_si,2)
     PHI(i,:) = PHI(i,:) + ...
                    p_si(1,j)*(meshg.pg(1,:).^p_si(2,j)) .* ...
                              (meshg.pg(2,:).^p_si(3,j));
   end
 end

 h = figure;
 pdemesh(h,grid,'k');
 hold on
 
 X = zeros(grid.ne,size(meshg.pg,2));
 Y = zeros(grid.ne,size(meshg.pg,2));
 UU = zeros(grid.ne,size(meshg.pg,2));
 
 for ie = 1:grid.ne
    
   XY_plot = grid.e(ie).bk * meshg.pg;
   XY_plot(1,:) = XY_plot(1,:) + grid.e(ie).xyv1(1);
   XY_plot(2,:) = XY_plot(2,:) + grid.e(ie).xyv1(2);

   U_el = NU((ie-1)*base.pk+1 : ie*base.pk);
   U_plot = U_el' * PHI;

   X(ie,:) = XY_plot(1,:);
   Y(ie,:) = XY_plot(2,:);
   UU(ie,:) = U_plot.*exp(xi(XY_plot)-xie(ie));
   
 end

 for ie = 1:grid.ne
   pdesurf_oct(h,[X(ie,:);Y(ie,:)],meshg.tg,UU(ie,:));
 end
 title('u_h^*');
 axis equal

 hold off
 
return

